Dark-bright solitons in Bose-Einstein condensates at 
finite temperatures 



V. Achilleosi, D. Yan^, P.G. Kevrekidis^, and D.J. 
Prantzeskakis^ 

^ Department of Physics, University of Athens, Panepistimiopohs, Zografos, Athens 
15784, Greece 

^ Department of Mathematics and Statistics, University of Massachusetts, Amherst, 
MA 01003-4515 USA 

Abstract. 

We study the dynamics of dark-bright sohtons in binary mixtures of Bose gases at 
finite temperature using a system of two coupled dissipative Gross-Pitaevskii equations. 
We develop a perturbation theory for the two-component system to derive an equation 
of motion for the soliton centers and identify different temperature-dependent damping 
regimes. We show that the effect of the bright ("filling") soliton component is to 
partially stabilize "bare" dark solitons against temperature-induced dissipation, thus 
providing longer lifetimes. We also study analytically thermal effects on dark-bright 
soliton "molecules" (i.e., two in- and out-of-phase dark-bright solitons), showing that 
they undergo expanding oscillations while interacting. Our analytical findings are in 
good agreement with results obtained via a Bogoliubov-de Gennes analysis and direct 
numerical simulations. 
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1. Introduction 

Macroscopic nonlinear excitations of atomic Bose-Einstein condensates (BECs) [HE] 
have been a subject of intense theoretical and experimental research over the last 
few years |3J. More specifically, matter- wave dark and bright solitons, that can be 
formed in single-component BECs with repulsive or attractive interatomic interactions 
respectively, have been observed in a series of experiments while their statics and 
dynamics have been extensively studied theoretically in various settings (see, e.g., [IHH] 
for recent reviews). Of particular interest are coupled dark-bright (DB) solitons that 
may exist in binary mixtures of BECs with repulsive interatomic interactions (such as 
ones composed by different hyperfine states of ^^Rb atoms [3 [8]): these solitons are 
frequently called symbiotic ones, as the bright soliton component (which does not exist 
in the system with repulsive interactions [1]) can be supported due to the nonlinear 
coupling with the dark soliton component. Such structures have recently been observed 
experimentally in a ^''Rb BEC mixture using a phase-imprinting method |9] or in 
two counter-flowing ^^Rb BECs [10l[TT], while they have also been studied in various 
theoretical works in continuum [TOHTS] and discrete [H] settings. 

The above theoretical studies on atomic DB solitons have been performed in the 
ideal case of zero temperature: in fact, finite-temperature induced dissipation of matter- 
wave solitons have basically been studied, so far, in the simpler case of dark solitons 
in single-component BECs [T5H20]. In particular, this problem was first addressed in 
Ref. [I5] (see also Ref. [IE]), where a kinetic-equation approach, together with a study 
of the Bogoliubov-de Gennes (BdG) equations, was used. In that work, it was found 
that the dark soliton center obeys an equation of motion of a harmonic oscillator, which 
incorporates an anti-damping term accounting for the finite temperature effect. The 
presence of this term alters the soliton trajectories so that the experimentally observed 
dark sohton dynamics can be qualitatively understood: solitons either decay fast at 
the rims of the BEC (for high temperatures) [211423] or perform oscillations of growing 
amplitude (for low temperature) [9l [24]426] and eventually decay. A similar equation 
of motion for the dark soliton center was also derived in Ref. [17] by applying the 
Hamiltonian approach of the perturbation theory for dark matter-wave solitons [6] to 
the so-called dissipative Gross-Pitaevskii equation (DGPE). This model incorporates a 
damping term (accounting for finite temperature), first introduced phenomenologically 
by Pitaevskii [27], and later shown to be relevant from a microscopic perspective (see, 
e.g., the review [28]). It is important to note that, as shown in Ref. [17], the analytical 
results obtained in the framework of the DGPE were found to be in very good agreement 
with numerical results obtained in the framework of the stochastic Gross-Pitaevskii 
equation (SGPE); see, e.g., Ref. [29] for a review on the SGPE model. It should also be 
mentioned that while the above works chiefly considered finite temperature effects for 
the case of a single dark soliton, the DGPE model and the anti-damping-incorporating 
ordinary differential equations (ODEs) for the soliton center were also examined in the 
case of multiple dark solitons. In particular, the cases of two and three oscillating and 
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interacting, anti-damped dark solitons were considered in Ref. |30] . 

In the present work, we study finite-temperature dynamics of DB solitons in 
liarmonically confined Bose gases. In particular, we adopt an effective mean- 
field description and analyze theoretically and numerically a system of two coupled 
DGPEs, describing the evolution of a binary quasi-one- dimensional (ID) BEG at finite 
temperature. We extend the considerations of Ref. [17] and develop a Hamiltonian 
perturbation theory for the two-component system at hand. This way, we obtain an 
equation of motion for the DB soliton center, similar to the one derived in Refs. [ISIITT] . 
This equation, which includes an anti-damping term accounting for finite temperature, 
provides a characteristic eigenvalue pair (i.e., a pair of solutions of the characteristic 
equation associated with the linear equation of motion), which is connected to the 
eigenvalue associated with the anomalous mode of the DB soliton. Performing a 
Bogoliubov-de Gennes (BdG) analysis, we show that the anomalous mode eigenvalue 
becomes complex as the dissipation (temperature-dependent) parameter is introduced, 
leading to an instability of the DB soliton pair. The temperature-dependence of the 
eigenvalues (determined analytically) is found to be in good agreement with the one of 
the anomalous mode eigenvalue (determined numerically). 

Furthermore, these considerations are generalized in the case of a DB soliton 
"molecule", composed by two-DB-solitons. In the latter setting, both configurations 
featuring in-phase and out-of-phase bright components can be obtained in the trap |31] . 
We illustrate their dynamical instabilities as a function of temperature and capture them 
analytically by means of coupled nonlinear DDEs accounting for the three ingredients 
(trap restoring force, interaction between DB solitons and thermally induced anti- 
damping). We show that, due to finite temperature, the nature of their interaction 
(and collisions) changes: for short times individual solitons behave as repelling particles, 
while for longer times they gain kinetic energy and completely overlap at the collision 
point. Our analytical considerations and numerical results reveal a fundamental effect: 
the partial stabilization that the bright ("filling") soliton component offers to the 
corresponding "bare" dark soliton against temperature-induced anti-damping. This way, 
a significantly longer lifetime of the symbiotic (dark-bright) structure can be achieved, 
in comparison to its bare dark soliton counterpart. 

The paper is structured as follows. In section II we present the model and study 
some of its basic properties such as the evolution to the equilibrium state. In section 
III we develop the perturbation theory to derive and solve the equation of motion for 
the single DB soliton; we also compare our analytical findings to numerical results. In 
section IV we generalize relevant considerations to the case of multiple DB solitons, and 
in section V we present our conclusions. 
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2. The model and its basic properties 



2.1. The system of dissipative Gross-Pitaevskii equations 

We consider a two-component elongated (along the x-direction) repulsive Bose gas, 
composed of two different hyperfine states of the same alkali isotope, and confined in a 
highly anisotropic trap (such that the longitudinal and transverse trapping frequencies 
are Ux <^ uj±). In such a case, the system can be considered as quasi-lD and, hence, 
the coupling constants take their effectively ID form, namely gjk = 2hw±ajk, where 
ttjk denote the three s-wave scattering lengths (note that ai2 = a.21) which account for 
collisions between atoms belonging to the same (ajj) or different {ajk,j 7^ k) species. Let 
us now focus on the experimentally relevant case of a two-component BEC consisting 
of two different hyperfine states of ^^Rb, such as the states 1) and |2, 1) used in 
the experiment of Ref. [8], or the states |1, —1) and |2, —2) used in the experiments of 
Refs. [inilll]. In the first case, the scattering lengths take the values an = 100. 4ao, 
ai2 = 97.66ao and 022 = 95.00ao, while in the second case the respective values are 
an = 100. 4ao, ai2 = 98.98ao and a22 = 98.98ao (where is the Bohr radius). In either 
case, it is clear that the scattering lengths and, accordingly, the effectively ID coupling 
constants take approximately the same values, say a^j ~ a and gij ^ g = 2hiJ±a, 
respectively, which is what we will assume henceforth. 

We now consider the case where the two-component Bose gas under consideration 
is at finite temperature. In particular, we assume that the thermal modes of energies 
> hw± are at equilibrium, accounting for a heat bath in contact with the axial part 
of the gas, while the modes in the x-direction are highly occupied so that the classical 
field approximation is valid [32','33] . Then, extending considerations pertinent to single- 
component Bose gases [281129 .32 .33 ] to the two-component case, we may use the following 
set of two coupled ID SGPEs to describe the axial modes of the system: 

ihdt^j = [1 - 7,(x, t)] (^-^dltPj + vix) -f^j+gYl l^^l' j + ^j-^^' 

Here, ipj{x,t) {j = 1,2) are complex order parameters characterizing each component 
of the binary Bose gas, m is the atomic mass, fj,j are the chemical potentials, 
while V{x) = {l/2)mujlx'^ is the external trapping potential. Furthermore, rij{z,t) 
are complex Gaussian noise terms with correlations of the form {ri*[x,t)r]j{x' ,t')) = 
2h'yj{x, t)kBT6{x—x')6(t—t'), where brackets denote averaging over different realizations 
of the noise. The strength of the latter can be calculated ab initio by the Keldysh self- 
energy [32]; for thermal clouds close to equilibrium, the relevant integrals determining 
the dissipation 7^(0;, t) can be expressed as follows: 

X [iVi(l + JVaXl + N,) + (1 + Ni)N,_N,], (2) 

where /3 = l/ksT, are the condensate energies, eiP are the energies of the n-th 
excited states, M^'^ = [exp(/3(Ei-''^ + 1/(a;)+2^X;Li(lV'fcP)-Aij))-l]"^ are Bose-Einstein 
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distributions, while En and /c„ = ^/2mEn/h? denote, respectively, the kinetic energies 
and momenta of single particles in the n-th excited state. Physically speaking, Eq. ([2]) 
describes the exchange of atoms between the thermal clouds and the condensates due 
to elastic collisions; notice that in the above description we have taken into regard 
exchanges up to the third excited state while, to leading order approximation, we 
have omitted exchanges between the different hyperfine states (in other words, we have 
considered the simplest situation where each condensate component interacts with its 
own thermal cloud). 

Under the above assumptions, the dissipation terms 7j(a;) may in principle be 
calculated numerically, for several temperatures, as was done in the case of a single- 
component Bose gas in Refs. |17]. In this work, it was shown that, sufficiently close 
to the trap center (i.e., in the interval [—R/2, R/2], where R is the Thomas-Fermi 
radius), the dissipation takes approximately constant values for a relatively wide range 
of temperatures. Furthermore, as shown in Refs. [15l[16l[T9] (see also the discussion in 
Ref. [I71[T9] and, more recently, in Ref. [20]), the value of 7 — which determines the dark 
soliton's life time — scales with temperature as 7 oc T"", with 1 < a < 4; note that the 
case 7 oc T'* corresponds to the regime ksT ^ /i, while the case 7 oc T corresponds to 
the regime fc^T ^ /i (where /i is the chemical potential of the background Bose liquid). 

Taking into regard the above findings, below we will consider the situation where 
both dissipative terms jj are constant: such an assumption is consistent with our scope, 
i.e., to analyze the dynamics of the DB-soliton near the center of the trap. Furthermore, 
based on the fact that simulations investigating soliton dynamics in the framework of the 
SGPE model were found to be in fairly good agreement with analytical and numerical 
results relying on the respective DGPE model, below we will omit the noise terms 
rij{x,t); this way, we will use the following system of two coupled DGPEs to describe 
the DB soliton dynamics in the two-component Bose gas at finite temperatures: 

{z - ^,)hdtij, = (^-^dli:, + V{x) -^^^+gJ2 l^fcl' j (3) 

Note that the above model was recently used in Ref. p4], where the quantum Kelvin- 
Helmholtz instability of a two-component BEG was studied. 

The system of Eqs. (|3]) can be expressed in dimensionless form as follows. Measuring 
the densities iV'jP, length, time and energy in units of 2a, a± = \ffifjJ^, u^^^ and hxi^^ 
respectively, Eqs. ([3]) become: 

{i - jd)dtUd = - ^dlud + V{x)ud + {\ud\^ + IwfoP - fi)ud, (4) 

{i - lb)dtUb = - ^dlub + V{x)ub + {\ub\^ + - /i - A)ub, (5) 

where we have used the notation ipi = Ud and il)2 = Ub, indicating that the component 1 
(2) is supposed to support a dark (bright) soliton, and the respective chemical potentials 
are now ni = Hd = H and = /it = /i + A; in our considerations below we assume that 
Ud > fJ'b, i-e., A = — |A| < 0. Finally, the external potential in Eqs. ©-(IH]) takes the 
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form V{x) = (l/2)i7^x^, where Q = Ux/u)± ^ 1 is the normahzed trap strength; the 
latter, along with the thermally induced damping parameters jdfi, are considered to be 
small parameters of the system (these will be treated as formal perturbation parameters 
in our analytical approximation - see below). 

We should add a comment here about the relevant range of values of the 
parameter 7. A number of recent experiments, including the ones in Hamburg [9||25]. 
Heidelberg pHl26] and Pullman p^jlTT] , have focused on regimes of very low temperature 
where the effect of the term associated with 7 is imperceptible (over the experimentally 
relevant time scales). The focus of these experiments was on the soliton dynamics 
and an effort was made (by operating at T/Tc < 0.1) to correspondingly minimize 
the thermal effects. It is easier to appreciate the latter features in the context of the 
earlier experiments of the Hannover group pTll23]. which were conducted in the regime 
of T/Tc ~ 0.5. In that realm, the relevant values of 7 can be estimated to be up to 
10"^ ^35j. In what follows, we will treat 7 generally as a free parameter, in order to 
illustrate the available wealth of bifurcation and dynamical phenomena of this system. 
Nevertheless, the reader more keen on the physical applications of the model to the 
physics of finite-temperature BECs should keep in mind the above values as a guideline 
towards the parameter regimes pertinent therein. We finally note that our analysis may 
also be used as a theoretical basis for understanding results of future experiments on 
dark and dark-bright solitons exploring finite-temperature effects (see, e.g., discussion 
in the Supplemental Material of Ref. |10j). 




2.2. Relaxation to the ground state of the system 

Since our purpose is to study the dissipative dynamics of DB solitons in this setting, it 
is natural to consider at first the dynamics of the pertinent background wave functions, 
namely a Thomas-Fermi (TF) wave function for the Ud component and a zero wave 
function for the Ui, component. In particular, we will show that the coupled DGPEs 
Eqs. (!4])-(l5]), similarly to their one-component counterpart (see, e.g., discussion in 
Ref. [36]), describe a relaxation process. Namely, as a result of the finite temperature, 
the two components, starting (at t = 0) from suitable initial conditions, will evolve so 
that, at sufficiently large times, Ud will converge towards a TF cloud with the prescribed 
value of the chemical potential while Ub will vanish. 

To show that this is the case indeed, we examine the peak amplitudes Ud,b(t) of 
the wave functions Ud^i^ = 0,t), corresponding to their (absolute) values at the center 
of the trap (i.e., at x = 0, where V{x) = as well), and assume respective phases 
6dfi{t). The evolution equations for Ud^it) and Od^bit), which can directly be obtained 
by introducing the ansatz Ud^ = Ud^it) exp[— 'i6'(i^b(t)] into Eqs. (|l])-([5]), are of the form: 




Ud,h + ld,bUd,b6d,b = 0, 

7dUd - 0dUd + {Ul + Ul - fi)Ud = 

IbUb - ObUb + {Ul + U^-fi-A)Ub = 



(6) 
(7) 
(8) 
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where overdots denote time derivatives. Next, utilizing Eqs. we obtain from Eqs. ([7])- 
([8]) the following system: 

tfd= - Id {U! + Ul - /i) f/rf, (9) 
f/fc = - 7b [Ul + U^-fi-A) U,, (10) 

where '^d,b = 7d,b/(l + lib)- It is clear that that the system of Eqs. (!9l)- (fT0l) has a 
fixed point {Udo, U^q) = {^/Jl, 0) [a similar analysis can be done for the fixed point 
{Udo, Ubo) = (0, v^/i + A)]. The evolution of small perturbations Udi,bi around this 
fixed point can then readily be found introducing the ansatz Udo{t) = ^/Jl + Udiif) 
and f/foo = Uhi(t) into Eqs. (l9ll-(fT0l) and linearizing with respect to Udi,bi] this way, we 
can easily solve the equations for Udi^i and finally obtain the following approximate 
expressions for the peak amplitudes of the wave functions: 

Udit) {Ud{0) - (11) 

Ub{t) ^ t/,(0)e-^^l^l*, (12) 

where Ud^biS^) are initial conditions. Thus, at sufficiently large times, the peak amplitude 
of Ud will decay to the value ^JJi, while the one of Ub will become zero. Accordingly, 
during the relaxation to equilibrium process, one may expect the following type of 
evolution towards relaxation. If the Ud component is initially a Thomas-Fermi (TF) 
cloud of amplitude ?7d(0), its density will evolve as, 

\ud{x,t)\^^Ul{t)-V{x). (13) 

On the other hand, if the Ub component has initially the form of an arbitrary localized 
function, e.g., a Gaussian, of amplitude f/b(0), it will asymptotically approach the trivial 
stationary state. 

The above predictions can be directly compared to numerical simulations. In 
particular, in Fig.[T]we show the evolution of a state characterized by the initial densities 
|Mrf(a;,0)|2 = Uj{Qi) - (l/2)nV and |Mfe(x,0)p = Ul{Q)eyi^[-2{x/d)\ with parameter 
values ?7d(0) = 0.86, f/fe(0) = 0.6, VL = 0.05 and d = 10; as found by direct numerical 
integration of Eqs. (|l])-([5]), with fi = 1.3, |A| = 0.1 and ■jd = lb = 0.05. The figure 
clearly shows the validity of our analytical approximations: the Ud component develops 
into a TF cloud with chemical potential n = 1.3, with the numerically found density 
profile [solid (red) line] being in fairly good agreement with the analytical prediction of 
Eq. ( lT3l) (dashed line); on the other hand, Wfc-component [solid (green) fine] vanishes at 
if: ^ 200, a time consistent with the slow time scale = (7f,|A|)~-^ ^ 200 suggested by 
Eq. dn]). 



3. Dissipative Dynamics of a Single Dark-Bright Soliton 

3.1. Analytical results 

Having studied the relaxation process described by Eqs. (|4])-(l5]), we will now proceed 
to investigate, in the same framework, the dissipative dynamics of DB solitons. We 
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Figure 1. (Color online) Time evolution of a state characterized by the densities 
Udix,0)\^ = C/J(0) - (l/2)l]2a;2 and \ub{x,0)\^ = [/^^(O) exp [-2(j;/w)2] , with 
parameter values Ud{0) = 0.86, Ub{0) — 0.6, — 0.05 and w = 10. The solid lines 
show the density of the dark (red) and bright (green) component, while the dashed line 
shows the analytical result of Eq. (fTTj). The other parameter values used in Eqs. Q-® 
are /i — 1.3, |A| — 0.1 and jd = Ih = 0.05. 



will assume that the dark soliton is on top of an already formed TF cloud with the 
equilibrium density |Md,TFp = fi' — V{x); this way, the density \ud\'^ in Eqs. (!l])-(l5]) is 
substituted by \ud\^ |Md,TFp|wdp. Furthermore, we introduce the transformations 
t ^ fit, X ^ y/Jlx, \ub\'^ — 7- yU^^l^bp, and cast Eqs. (HJ-dS]) into the following form: 

idtUd + ^dlud - {\ud\^ + -l)ud = Rd, (14) 

idtUb + ^dlub - {\ub\^ + ludp - fi)ub = Rb, (15) 
where fl = 1 + A//i and 

Rd = (2^2)-i[2(l - \ud\^)Vix)ud + V'ix)d^Ud] + idfi-'dtUd, (16) 

Rb = Ai"^[(l - \ud\'^)V{x)ub + fi'^bdtUb]. (17) 

while V'{x) = dV/dx. Equations fll4p - (IT^ can be viewed as a system of two 
coupled perturbed nonlinear Schrodinger (NLS) equations, with perturbations given 
by Eqs. (IT6l)-( fT7|) . In the absence of the perturbations, i.e., at zero temperature 
(76 = 7d = 0) and for the homogeneous system iV{x) = 0) subject to the boundary 
conditions Iw^p — )■ 1 and — >■ as |a;| ^ oo, the NLS Eqs. (IT^- (fT5l) possess an exact 
analytical one-DB-soliton solution of the following form: 

Ud{x,t) = cos0tanh[D(x — xo(t))] + isin0, (18) 
Ub{x, t) = r7sech[D(x - xo(t))]exp[ikx + i6'(t)], (19) 

where is the dark soliton's phase angle, cos0 and t] represent the amplitudes of the 
dark and bright solitons, and D and Xo{t) are associated with the inverse width and the 
center position of the DB soliton. Furthermore, k = Z?tan0 = const and 6{t) are the 
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wavenumber and phase of the bright sohton, respectively. The above parameters of the 
DB-sohton are connected through the following equations: 

D'^ = cos^ (p - rf, (20) 
x'o = L'tanc/), (21) 

e{t) = \{D'-k^)t + {^/^,% (22) 

with xq denoting the DB soliton velocity. Notice that the amplitude i] of the bright 
soliton, the chemical potential /i of the dark soliton, as well as the (inverse) width 
parameter D of the DB soliton are connected to the number of atoms of the bright 
soliton by means of the following equation: 

N,^ ( \u,\'dx = (23) 

Let us now employ the Hamiltonian approach of the perturbation theory for the 
matter-wave solitons to study the dissipative dynamics of DB solitons. We start by 
considering the Hamiltonian (total energy) of the system of Eqs. ( IT^ - (JT5l) . in the absence 
of the perturbations (i.e., for Ri, = Rd = 0), namely. 
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Loo 



E = - £dx, 

£ = \d,Ud\^ + \d,Ub\^ + {\ud\^ + \ub\^ - I f - 2{fi - l)\ub\^. (24) 

The energy of the system, when calculated for the DB soliton solution of Eqs. f fT8|) -f fT9|) . 
takes the following form: 



3 V2 l^J 

We now consider an adiabatic evolution of the DB soliton and, particularly, we assume 
that, in the presence of the perturbations of Eqs. f|T6|) -f lT7|) . the DB soliton parameters 
become slowly-varying unknown functions of time t. Thus, the DB soliton parameters 
become — )■ 0(t), D — )■ D{t) and, as a result, Eqs. (!20|) -( l2T|) read: 

D^{t)= cos" <P{t)-^D{t), (26) 

Xo(t) = D(t)tan0(t), (27) 

where we have used Eq. f l2^ . The evolution of the parameters 0(t), D{t) and XQ(t) can 
be found by means of the evolution of the DB soliton energy. In particular, employing 
Eq. (12^ . it is readily found that 
dE 

— = 4DD^ + xD sec^ (j){D + D0tan0) . (28) 



On the other hand, using Eqs. (]T4l) - (]T5l) and their complex conjugates, it can be found 
that the evolution of the DB soliton energy, due to the presence of the perturbations, 
is given by: 

^ = - 2Re jy {R*Aud + RldtUb)dx'^ , (29) 
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where the asterisk denotes complex conjugation. Substituting and Rb into Eq. 
and evaluating the integrals, we finally obtain from Eqs. f l28|) -f l29|) the following result: 

4:DD^ + xD sec^ (j){D + L'c/.tan^) = 

— (2 cos^ sin — xD sin cos 0) V'(xo) 

_^l£D'sm^(j)--^xDHan'(l). (30) 
3 /i 3 /i 

Equation (!30|) . together with Eqs. ( |26|) - (!27|) . constitute a system of equations for the 
unknown soliton parameters 0(t), D{t) and xo(t). In the case of a DB soliton near 
the center of the trap with an almost "black" dark-soliton-component (i.e., Xq ~ and 
COS0 ~ 1), the above system has a fixed point Xo^eg = and (peg = 0, and 



Considering now small perturbations around the fixed points, i.e., xq — )■ O+xq, 4> — )■ + 
and D — > D^q + Di, we linearize Eqs. fl2Bl) - fl7n) and Eq. fl5Ul) with respect to Xq, and 
Di, and obtain the following results: 

-1 

2, 



D,= -b<p\ D= (2D,q + ^) \ (32) 



-2 + X^eg 



V'ixo) 



D,q[-SD,qb-x{2D-D,q)]' 

xo = Deg<p. (34) 

Differentiating Eq. ( I34l) with respect to time once, and using Eqs. (I33l)- (1M1) . we obtain 
after some straightforward algebraic manipulations the following equation of the motion 
for the DB soliton center xq: 

Xo - axo + w^sc^^o = 0, (35) 

where the oscillation frequency Wqsc and the anti-damping parameter a are respectively 
given by: 



<c = ^' Xo = 8Wl+^ , (36) 



1- 2 ^ 4 X / 1 2 



a = [1. - 8^ 7.J + 3-^ ^7. - 7. + gX 7.J . (37) 

We note that in the absence of dissipation [a = in Eq. ( 15^ ]. Eq. recovers the 
results of Ref. [12]: according to this work, if both components are confined in the same 
harmonic trap of strength VL then a DB-soliton oscillates around the trap center with 
the frequency cUosc, given in Eq. 
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It is clear that the nature of the sohton trajectories Xo(t) as predicted by Eq. (15^ 
depend on whether the roots of the auxihary equation s"^ — as + cu^g^. = are real or 
complex. The roots are given by 



with the discriminant V = — a^^, determining the type of the motion. In particular, we 
identify different temperature-dependent damping regimes: the subcritical weak anti- 
damping regime (P < 0, a < Ocr), the critical regime {V = 0, a = acr), and the 
super-critical strong antidamping regime ("D > 0, a > acr)- In the first regime the 
soliton performs oscillations of growing amplitude, with xo(t) oc exp{at) cos(u;osc^), while 
in the latter two regimes the soliton follows an exponentially growing trajectory, i.e., 
Xo(t) oc exp(si^2^) (with si^2 ^ '^), and decays at the rims of the condensate cloud (see 
also below). 

3.2. Numerical results 

We now turn to a numerical examination of the above findings. First, we will show 
that our analytical predictions are supported by a linear stability analysis around the 
stationary DB sohton, say Uq = (u, v)-^ [see Eqs. (!T8!)- (fT9l) for = and xq = 0]. For 
such a state, the right-hand side of Eqs. Q-dS]) still vanishes and, thus, stationary DB 
solitons are exact solutions of the problem with 7^ ^, 7^ 0. We obtain this solution by 
means of a fixed point algorithm and then find the linearization spectrum around the 
stationary DB soliton state as follows. We introduce the ansatz 



into the DGPEs Eqs. (H])-© (here {A, (a, b)} define an eigenvalue-eigenvector pair, and 
e is a formal small parameter), and then solve the ensuing BdG eigenvalue problem. 

In Fig. [2], we observe a prototypical realization of a stationary DB soliton in a trap 
of strength = 0.1 (for simplicity, we consider the case with 7d = 76 = 7)- Notice 
that upon the variations of 7 (and hence of temperature) considered in the figure, 
the solution profile does not change, as mentioned above; however, the linearization 
problem and its eigenvalues significantly depend on the value of 7, as is shown in the 
four bottom panels of Fig. [2l In the zero-temperature (Hamiltonian) case of 7 = 0, all 
eigenvalues are imaginary. Furthermore, the oscillatory motion of a single DB soliton 
in the trap [12] (see also recent work in Refs. [IOl|Tl],|3T]) is spectrally associated with 
the existence of a single anomalous (alias negative Krein sign, "translational" ) mode in 
the linearization around the stationary soliton. In analogy to the case of dark solitons 
(see e.g. Refs. [I5l[26l|30]), this anomalous mode possesses a frequency identical to the 
frequency of the DB soliton oscillation, i.e., cuam = Ihi(Aam) = i^osc- 

Our analytical approximation for Wqsc is tested against the numerical results for 
Im(AAM), both in the case examples of Fig. [21 as well as in the parametric dependence 
results of Fig. [31 It is clear from the spectral plots (middle and bottom panels of Fig. [2D 
that, as soon as 7 7^ 0, the relevant anomalous eigenmode (indicated by red stars 




(38) 



u = uo + e[exp(At)a(a;) + exp(A*t)b*(x)], 



(39) 
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Figure 2. (Color online) The top panel depicts the stationary solution for a single 
DB-soliton for fi = 1.5, |A| = 0.6 and = 0.1. The dark (bright) components are 
shown by the dashed green (solid blue) lines. The middle and bottom panels are four 
spectral planes, corresponding to different values of "f — jd — Jb, for the single dark- 
bright soliton stationary states: in the middle left panel 7 = (the zero-temperature 
Hamiltonian case), in the middle right 7 = 0.05, in the bottom left 7 = 0.12, and 
in the bottom right 7 — 0.17. The (red) stars highlight the anomalous mode (of the 
Hamiltonian case) eigenvalues. 



in Fig. [3]) becomes complex, leading to soliton oscillations of growing amplitude; this 
behavior, which corresponds to the "subcritical" regime mentioned above, is similar to 
the case of dark solitons [171 130] and in accordance with rigorous results pertaining 
to dissipative NLS systems [37]. If 7 is increased beyond a critical point, namely 
jcr ~ 0.141, the relevant eigenvalue pair collides with the real axis, leading to the 
emergence of a pair of real eigenvalues (cf. bottom right panel of Fig. |2] and Fig. |3]). 
This corresponds to the "super-critical" regime mentioned above (see also Refs. [T7t[30]). 
where the divergence of the soliton from its center equilibrium is purely exponential. 
Notice that the analytical predictions for the relevant unstable eigenvalue (and the 
oscillatory or purely exponential divergence from the equilibrium position) in Figs. [2] 
and |3] are generally fairly accurate, although their accuracy is decreasing as 7 gets 
larger; this can be understood by the fact that our analytical approximation relies on 
the smallness of 7 which was treated as a small parameter of the problem within our 
perturbation theory approach. 

Last but not least, the role of the bright-soliton component in the dynamics 
should be highlighted in connection to the case of a dark soliton in a single-component 
condensate (where the bright soliton is absent). It can be directly seen from Eq. (137]) 
that the anti-damping effect is always weaker for the DB soliton in comparison to the 
dark one (at least in the case = lb = 1 consider herein). Hence, the lifetime of the 
DB soliton is always longer than that of the dark soliton and, in fact, it becomes larger. 
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Figure 3. (Color online) The real part (top) - instability growth rate- and imaginary 
part (bottom) -oscillation frequency- of the unstable eigenmode of the linearization 
around a stationary DB soliton as a function of the parameter 7 (parameter values are 
as in Fig. [2]). Solid (blue) lines indicate the full numerical result, while dashed (red) 
ones the analytical results of Eq. ((38|) . For 7 < jcr ~ 0.141 (subcritical regime), the 
complex conjugate pair is responsible for soliton oscillations of growing amplitude. For 
7 ^ 7cr (critical/super-critical regimes) the collision of the complex conjugate pair of 
eigenvalues creates a real pair, and the dynamics involves purely exponential growth. 
The parameter values are fj, = 1.5, |A| = 0.6, and fl = 0.1. 




Figure 4. (Color online) Comparison of the bifurcation diagrams Xr{j) for a 
stationary DB soliton [solid (blue) line] and a "bare" dark soliton [dashed (black) 
line] as obtained by the BdG analysis. It is clear that the presence of the bright 
("filling") component drifts 7cr towards larger values, acquiring also smaller values of 
the instability growth rate as compared to the ones found for the dark soliton. The 
parameter values are ft — 0.1, /i = 1.5 and |A| = 0.1. 
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Figure 5. (Color online) Contour plot showing the space-time evolution of the density 
in the single DB soliton case. The top panel represents the dark soliton and the bottom 
one the bright soliton, with 7 = 0.02 (subcritical regime) and ft = 0.1. The soliton 
is initially placed at a::o(0) = 0.4. The dashed line represents the analytical result 
of Eq. (1551) . namely xo{t) = a;o(0) exp[(a/2)t] cos(wosci). The parameter values are 
fi = 1.5, |A| = 0.6, and n = 0.1. 



as the bright-soliton component "filling" of the dark one becomes stronger. This partial 
stabilization of the dark soliton evolution by means of its symbiotic second component 
is clearly illustrated in Fig. |U where the bifurcation diagrams for the DB soliton are 
directly compared to the ones corresponding to the "bare" dark soliton. It is clear that 
the whole bifurcation diagram for the DB soliton is "drifted" towards larger values of 
7 (e.g., 7cr = 0.212 for the DB soliton and 7cr = 0.155 for the dark soliton), acquiring 
also smaller values of the instability growth rate A,- as compared to the ones found for 
the dark soliton for the same values of the temperature parameter 7. This is a clear 
indication that the "filled" dark soliton in a two-component BEC is more robust in the 
presence of finite temperature than a "bare" dark soliton in a single- component BEC. 
This is one of the principal findings of the present work. 

Our analytical predictions were also tested against direct numerical simulations 
illustrating the evolution of the single DB soliton, both for the sub-critical case of 
oscillatory growth (see Fig. |5]), and for the supercritical case of purely exponential 
growth (see Fig. |6]). In both cases it can be seen that the dashed line corresponding to 
the analytical solution of the ODE (!35|l accurately tracks the evolution of the center of 
the DB soliton, which progressively loses its contrast and eventually disappears in the 
condensate background, with the system converging to its ground state (see section 2). 

It is worth noting that in the results of Figs. |5] and |6] we have used, as initial 
condition, a TF cloud with a density at the trap center equal to the chemical potential 
fi appearing in the DGPEs (I4])-(I5]). Nevertheless, we have also briefly studied a case 
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Figure 6. (Color online) Similar to Fig. [5j but now for the super-critical case of 
7 = 0.15, which shows the exponential divergence of the soliton center. The soliton 
is initially placed at a;o(0) = 1. The dashed line represents the analytical result of 
Eq. namely xo{t) = [s2exp(sit) — siexp(s2t)]. 

where the density at x = of the TF cloud was different from /i. The evolution in such 
a far from equilibrium scenario, is shown in Fig. [TJ where the parameters are as in the 
subcritical case of Fig. |5l but with Ud{0) = 0.8. It is readily observed that apart from 
the transient period towards equilibrium (i.e., when the density is rearranged so that 
it properly corresponds to the relevant value of the chemical potential fi = 1.5), the 
agreement between analytical and numerical results is fairly good. I.e., the fast scale of 
the background relaxation does not substantially affect the evolution of the DB wave 
on the slower time scale of the oscillatory decay of the latter. 

4. Two Dark-Bright Soliton States 

We now focus on the study of DB soliton "molecules" composed by two-DB-soliton 
states. Let us first consider the homogeneous case {Q = 0), and use the following 
ansatz to describe a two-DB-soliton state composed by a pair of two equal-amplitude, 
oppositely located (at x = ±xo) single DB solitons: 

ipi{x,t) = (cos tanh X_ + i sin 0) (cos tanh X+ — i sin 0) , (40) 
Mx, t)=v sechX_ e'['=-+^(*)] + r] sechX+ e'I-^-+^(*)l e'^^ (41) 

where X± = D {x ± Xo{t)) , 2xo is the relative distance between the two solitons, and 
is the relative phase between the two bright solitons (assumed to be constant); below 
we will consider both the out-of-phase case, with A9 = vr, as well as the in-phase case, 
corresponding to = 0. Note that, similarly to the case of a single-DB soliton, the 
number of atoms A*";, of the bright-soliton component in the above two-DB-soliton state 
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Figure 7. (Color online) Same as the bottom panel of Fig. [5j but with an initial 
density of the TF cloud at the trap center equal to Ud{0) = 0.8. 



can be used to connect the DB-soliton parameters; in particular, if the two DB sohtons 
are well-separated then A^^^ is approximately twice as large compared to the result of 
Eq. ([2SD, namely, Nh ^ A'q'^^/D. 

As was recently shown in Ref. [31], at zero temperature (i.e., 7d = 76 = 0), the 
evolution equation for the DB soliton center (for f2 = 0) reads: 

(42) 



-^int, 

-^DD + -^BB + 2Fdb- 



(43) 



In the above equations, Fint is the interaction force between the two DB solitons, which 
consists of three different components: the interaction forces Fdd and Fbb between the 
two dark and two bright components, respectively, as well as the interaction force F^-q 
of the dark soliton of the one soliton pair with the bright soliton of the other pair (and 
vice- versa). These forces depend on the soliton coordinate Xq, as well as on the DB 
soliton parameters, as follows |3T] : 



DD 



BB 



+ 



1 

Xo 

x_ 



Xo 



-(544 - 352^2^) + 128De, {D% - l) xo 

- 6Deg + ^D%xo - 2x) D% cos A^e-^^'^'"" 
(1 + 2 cos^ A^) {-SDeqXQ + 6) 



4DegX0 



(44) 



Jj2 ^-ADeqXo 



eq 



(45) 

Fdb = ^(sz^egcos A^)e-2^-^o _ _ 64D,,xo)/^e,e-^''-"", (46) 

XO ^ ^ Yn V 3 / 



Xo 



where we have assumed that Z)(t) ~ and, thus, D{t) ~ D^q. 
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Figure 8. (Color online) The top panel depicts the stationary solution for an in- 
phase two-DB-soliton state for /i = 1.5, |A| = 0.6 and fl = 0.1. The dark (bright) 
components are shown by the dashed green (solid blue) lines. The middle and bottom 
panels are four spectral planes, corresponding to different values of 7, for the two DB- 
solitons in an in-phase, stationary configuration: in the middle left panel 7 = (the 
zero-temperature Hamiltonian case), in the middle right 7 = 0.1, in the bottom left 
7 = 0.2, and in the bottom right 7 = 0.4. 



Next, let us consider the case of two DB-solitons in the presence of the harmonic 
trap. Then, each of the two sohtons is subject to two forces: (a) the restoring force of the 
trap, Ftr [in the case of a single DB-soliton, this force induces an in-trap oscillation with 
a frequency cuosc — see Eq. ( |36l) ]. and (b) the pairwise interaction force Fint [cf. Eq. ( H3l) ] 
with other dark-bright solitons. Thus, taking into regard that Ft 



tr 



j2 „ 
osc 



write the effective equation of motion for the center xq of a two-DB-soliton state as 
follows: 

xo = Ftr + Fi„t. (47) 

In order to complete the consideration of the case at hand, we will finally study 
the finite-temperature effect on a two DB-soliton state in the trap. To do so, we will 
combine the thermal effect on each DB soliton in the trap, represented by Eq. ( l35l) . 
and interaction effects included in Eq. (142|) . This way, we may use the following 
approximation to describe the motion of the centers of the two DB solitons: 

Xo - axo - (Ft, + Fint) = 0. (48) 

The equilibrium points easily be found as solutions of the transcendental 

equation resulting from Eq. (148!) letting xq = Xq = in both the in- and out-of-phase 
cases. To study the stability of these equilibrium points in the framework of Eq. (HHj) . we 
use the ansatz xo{t) = Xgq + 6{t), and obtain a linear equation for the small- amplitude 
perturbation 6{t), namely: 6 — a6 + uf6 = 0, where the frequency Ui is given by. 



22,2 2 t^-^int 
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Figure 9. (Color online) The real (top panel) and imaginary (bottom panel) part of 
the two anomalous mode eigenvalues for a two in-phase DB soliton state. Both modes 
have complex, for 7 7^ 0, and both pairs eventually collide and give rise to exponential 
instabilities through eigenvalues on the real axis. Solid (blue) lines yield the numerical 
results, while the dashed (red) lines provide the corresponding theoretical predictions 
of Eqs. (p|-(j49|). 

where u^^^ and a are respectively given by Eq. (1361) and Eq. (jSTj). 

We now test the relevant predictions against BdG simulations, first for the in- 
phase case in Figs. |8]|9] and then for the out-of-phase case in Figs. ITOlfTTl As expected, 
in the case of the in-phase configuration the BdG analysis reveals the existence of two 
anomalous modes: the one with the smaller (larger) eigenvalue — in the zero-temperature 
case — corresponds to an in- (out-of-) phase motion of the two DB solitons, similarly to 
the case of a two-dark-soliton state in a single-component BEG [26|l30]. These two 
anomalous mode pairs lead to complex eigenfrequencies for 7 7^ 0, and the two-DB- 
soliton state performs oscillations of growing amplitude. Similarly to the case of the 
single DB soliton, as 7 is increased these pairs collide pairwise on the real axis in two 
critical points, namely 71 ~ 0.153 and 72 ~ 0.38. Beyond the second critical point 72, 
the growth of the trajectory of the DB soliton center becomes purely exponential. The 
theoretical approximation of the relevant complex (and subsequently real) eigenvalues 
depicted by dashed line in Fig. |9] is again fairly accurate, becoming progressively worse 
as 7 increases. 

A similar phenomenology arises in the case of out-of-phase two-DB-soliton states, 
as shown in Figs. [TOlfTTl However, there exists a rather nontrivial twist in comparison to 
the previous case. In particular, a third pair of complex eigenvalues emerges due to the 
fact that a third anomalous mode exists for 7 = 0. This mode is no longer a tmnslational 
one associated with the in- or out-of-phase motion of the two soliton centers (as before 
and as shown in the bottom left and bottom right eigenmodes of Fig. [T0|) . It is instead 
a mode associated with the tt relative phase of the peaks: if we add the eigenvector of 
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Figure 10. (Color online) The top panel depicts the stationary solution for two DB- 
solitons in an out-of-phase configuration, with parameters fi — 1.5, |A| = 0.6 and 
= 0.1. The dark (bright) soliton components are shown by the dashed green (solid 
blue) lines. The middle panels show three spectral planes, corresponding to different 
values of 7, namely from left to right we have 7 = (the zero-temperature Hamiltonian 
case), 7 = 0.02 and 7 = 0.2 respectively. In the bottom panel, we compare the bright 
soliton component's stationary solution (solid blue line), against the perturbed states 
(dashed red line) obtained by adding to it the respective BdG eigenfunctions. More 
specifically "A" , "B" and "C" correspond to the eigenfunctions of the three anomalous 
modes' eigenvalues in ascending order (see text). 

this unstable (for 7 7^ 0) mode to the two-DB-sohton out-of-phase solution, we observe 
that while the center location of the state remains intact, the relative heights of the two 
solitons are affected, leading to a symmetry breaking of the configuration. We will not 
consider this unstable mode further since its induced instability is weaker than those of 
the (in-phase and out-of-phase) translations. Nevertheless, we note that all three pairs 
of modes eventually collide on the real axis, eventually leading to pairs of purely real 
eigenvalues. 

Finally, we turn to direct numerical simulations for both the in-phase two-DB- 
soliton state in Fig. [121 and for the out-of-phase two-DB state in Fig. [131 In both cases, 
we show only the I0W-7, oscillatory growth (subcritical) regime. Despite the complexity 
of the resulting system and of the DB soliton interactions, it can still be clearly observed 
that the ODE (HHl) can be used to capture fairly accurately the relevant dynamics even 
for the long time evolutions considered in these figures. Here, it should be mentioned 
that the temperature- induced dissipation results in an interesting effect. Particularly, as 
observed in Figs. [121 and [131 fo^' short times, the individual DB solitons clearly behave 
like repelling particles, which can always be characterized by two individual density 
minima — even at the collision point. Nevertheless, for longer times, the nature of their 
interaction changes: due to dissipation, they gain kinetic energy and completely overlap 
at the collision point. 
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Figure 11. (Color online) The real part (top panel) and imaginary (bottom panel) 
part of the two anomalous mode eigenvalues (corresponding to "A" and "C" in the 
middle panel of Fig [TU)) for an out of phase two DB soliton state. Both modes lead 
to Hopf bifurcations, for 7 7^ 0, and both pairs eventually collide and give rise to 
exponential instabilities through eigenvalues on the real axis. The solid (blue) lines 
depict numerical results, while the dashed lines provide the corresponding theoretical 
predictions of Eqs. p5)) -(|i ^ . 




Figure 12. (Color online) Contour plot showing the space-time evolution of the 
density in the two dark-bright soliton in-phase state. The top panel represents the 
dark solitons and the bottom one the bright solitons, with 7 = 0.01. The solitons are 
initially placed at xi — 2.75 and X2 = —2.75. The dashed line represents the result 
obtained by numerical solution of Eq. pS)) . 
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Figure 13. (Color online) Contour plot showing the space-time evolution of the 
density in the two dark-bright soliton out-of-phase state. The representation of the 
solitons and the value of 7 is the same as in Fig. [121 The solitons are initially placed at 
xi = 1.76 and X2 — —1.76. The dashed line represents the result obtained by numerical 
solution of Eq. (|48|). 

5. Conclusions 

In the present work, we presented a systematic analysis of a prototypical model (the 
so-called dissipative Gross-Pitaevskii equation) incorporating the effects of temperature 
on the dynamics of dark-bright (DB) solitons. This was done both in the case of a single 
DB soliton, as well as in the case of DB soliton "molecules" , composed by multiple (in- 
or out-of-phase) DB solitons. 

We have developed a perturbation theory for the two-component system to 
analytically show the following: similarly to dark solitons, dark-bright ones execute 
anti-damped oscillations of growing amplitude for sufficiently low temperatures, while 
if the relevant parameter becomes sufficiently large, then the decay of the contrast of 
the solitons (and their disappearance in the background) becomes exponential. 

A fundamental effect revealed by our analysis is that the presence of the bright 
("filling") component hinders the temperature-induced dissipation associated with the 
dark soliton, and offers a significant partial stabilization (i.e., a significantly longer 
lifetime) to the corresponding symbiotic DB soliton structure, in comparison to its 
"bare" dark soliton counterpart. The above effect relies on the fact that the critical 
value of the relevant parameter (labeling the different damping regimes) is increased, 
while the instability growth rate is decreased, for the DB solitons. Similar conclusions 
were reached in the case of two dark-bright entities, with the added twist that their 
relative phase may introduce (in the out-of-phase case) additional anomalous modes 
and instability sources in the system. The latter are not associated with in- or out-of- 
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phase translational motion of the sohtons but rather with a symmetry-breaking in their 
relative amphtudes. 

As concerns the relevance of our findings with pertinent experimental efforts 
we note the following. First, all relevant recent experiments for dark and dark- 
bright solitons were conducted at extremely low temperatures, aiming to minimize 
corresponding anti-damping effects. In our setting this corresponds to subcritical 
dynamics of small 7. Nevertheless, we believe that our findings may be relevant to 
future experiments exploring in more detail finite-temperature effects (see Supplemental 
Material of Ref. [lO]). 

A natural direction to extend the present studies is to consider the higher 
dimensional setting of vortices [38] and of their two-component generalizations, namely 
the vortex-bright solitons |39]. Understanding the thermally induced dynamics and the 
modifications of the corresponding precessional motion, especially in the presence of 
multiple coherent structures would constitute an interesting topic for future study. On 
the other hand, it would certainly be relevant to extend the present studies to more 
complex models that provide coupled dynamical equations for the condensate and the 
thermal cloud [28] (rather than use a single equation directly incorporating the effects 
of the thermal cloud on the condensate without the possibility of "feedback"). Such 
studies are in progress and pertinent results will be reported in future publications. 
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